This notebook is Part 3 of an exploration to visualize results of hydrologic models. In Part 1 we built custom HTML widgets using D3.js, and in Part 2 we looked at rendering water fluxes using Sankey diagrams. Here we test multiple libraries to generate hillshade (3D) views of river basins and water infrastructure, in particular we want to compare leaflet, Three.js and MapboxGL implementations.
Aside from rendering spatial terrain and water streams in 3D (and potentially other covariate layers), our objective is to overlay custom labels to illustrate and quantify the hydrological cycle.
Some inspiration below:
- tylor Morgan Tutorial: Adding OSM Data to Rayshader Maps
- NASA Web WorldWind
- Flood depth on 3D terrain
- Contouring with Threes.js
Data Acquisition
We’ll start with a sample scene of the Selingue Dam in the Niger River basin.
basin <- shapefile("~/Projects/2021-iwmi/data/mli/srtm/mli_basin.shp")
zoi <- shapefile("~/Projects/2021-iwmi/data/mli/srtm/zoi.shp")
ext <- extent(zoi)
center <- coordinates(zoi)
# Get CGIAR SRTM DEM at 90m
getData("SRTM", lon=center[,1], lat=center[,2], path="_data") %>%
crop(zoi) %>%
mask(zoi) -> srtm
srtm
class : RasterLayer
dimensions : 838, 973, 815374 (nrow, ncol, ncell)
resolution : 0.0008333333, 0.0008333333 (x, y)
extent : -8.6725, -7.861667, 11.06917, 11.7675 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : srtm_35_10
values : 326, 476 (min, max)# Get LUC
#luc <- raster("~/Projects/2021-iwmi/data/mli/stat/")
# Admin boundaries
The basin covers a large area, so we need 8 SRTM tiles, but 1 is enough for a proof of concept. Next we’ll get a satellite basemap.
# Satellite basemaps
maptiles::get_tiles(terra::ext(zoi), "Esri.WorldImagery", zoom=10) %>%
stack() %>%
crop(zoi) %>%
mask(zoi) -> bmap
plot(terra::vect(basin), axes=T, asp=1, col=pal[2], border=pal[1], lwd=2,
main="Niger River Basin (Mali)")
plot(zoi, add=T, lty=3, col=alpha(pal["red"], .6), border=pal["red"], lwd=2)
text(-8, 10.5, "Selingue Dam\n(Mali)", col=pal["red"], cex=.7, font=2)
grid()
plot(terra::ext(zoi), axes=T, asp=1, border=NA,
main="Selingue Dam (Mali) - ESRI World Imagery")
plotRGB(bmap, add=T)
grid(col="white")
plot(srtm, axes=T, asp=1, border=NA,
main="Selingue Dam (Mali) - SRTM 90m")
grid(col="white")
WebGL (rayshader)
Next we convert the 2 rasters to a matrix format that’s compatible with Rayshader methods.
Show code
# Convert rasters to rayshader matrix format
srtm_array <- raster_to_matrix(srtm)
# Convert sat basemap to matrix (test)
r <- raster_to_matrix(bmap$red)
g <- raster_to_matrix(bmap$green)
b <- raster_to_matrix(bmap$blue)
bmap_array <- array(0, dim=c(nrow(r), ncol(r), 3))
bmap_array[,,1] <- r/255
bmap_array[,,2] <- g/255
bmap_array[,,3] <- b/255
bmap_array %>%
aperm(c(2,1,3)) %>%
# Stretch contrast
rescale(to=c(0,1)) -> bmap_array
srtm_water <- srtm_array
srtm_water[srtm_water < 353] <- 0
basemap_sat <- srtm_array %>%
height_shade() %>%
add_overlay(bmap_array) %>%
add_shadow(ray_shade(srtm_array, zscale=90)) %>%
add_water(detect_water(srtm_water), color=alpha(pal["blue"], 0.4))
basemap_sat %>% plot_map()
Figure 1: Hillshade Basemap of Selingue Dam (Niger River basin)
alt <- getData("alt", country="MLI", path="_data") %>%
crop(zoi) %>%
raster_to_matrix()
basemap_sat %>% plot_3d(alt,
zscale=20, zoom=0.5,
phi=30, theta=20, fov=0,
shadow=TRUE, shadowcolor=pal["black"])
rglwidget()
Figure 2: Hillshade Basemap of Selingue Dam (Niger River basin)
That doesn’t look very clear, so instead we’ll create a basemap, not using the satellite image but a simple desert texture.
base_map <- srtm_array %>%
height_shade() %>%
add_overlay(sphere_shade(srtm_array, texture="desert",
zscale=4, colorintensity=5), alphalayer=0.5) %>%
add_shadow(lamb_shade(srtm_array, zscale=6), 0) %>%
add_shadow(ambient_shade(srtm_array), 0) %>%
add_shadow(texture_shade(srtm_array, detail=8/10, contrast=9, brightness=11), 0.1) %>%
add_water(detect_water(srtm_water), color=alpha(pal["blue"], 0.4))
saveRDS(base_map, "./_data/base_map.rds")
Show code
Figure 3: Step 1 - Hillshade Basemap of Selingue Dam
Looks better, so let’s acquire and overlay spatial features from OSM.
osm_roads <- opq(bbox(zoi)) %>%
add_osm_feature("highway") %>%
osmdata_sf()
osm_water = opq(bbox(zoi)) %>%
add_osm_feature("waterway") %>%
osmdata_sf()
osm_place = opq(bbox(zoi)) %>%
add_osm_feature("place") %>%
osmdata_sf()
road_layer <- generate_line_overlay(
dplyr::filter(osm_roads$osm_lines, highway %in% c("primary", "secondary")),
extent=ext, srtm_array, linewidth=5, color=pal["black"])
water_layer <- generate_line_overlay(
osm_water$osm_lines,
extent=ext, srtm_array, linewidth=3, color=pal["blue"])
place_layer <- generate_label_overlay(
dplyr::filter(osm_place$osm_points, !is.na(name) & nchar(name)<10),
extent=ext, heightmap=srtm_array,
font=2, text_size=1.6, point_size=1.6, color=pal["black"],
halo_color="white", halo_expand=2, halo_blur=1, halo_alpha=.9, seed=1,
data_label_column="name")
scene <- base_map %>%
#scene <- basemap_sat %>%
add_overlay(road_layer) %>%
add_overlay(water_layer, alphalayer=1) %>%
add_overlay(place_layer)
saveRDS(scene, "./_data/scene.rds")
Finally we’ll use WebGL to render this scene in 3D.
amb_layer <- ambient_shade(srtm_array, zscale=1/5)
scene2 <- srtm_array %>%
height_shade() %>%
add_shadow(texture_shade(srtm_array, detail=8/10, contrast=9, brightness=11), 0) %>%
add_shadow(amb_layer, 0) %>%
add_overlay(road_layer) %>%
add_overlay(water_layer, alphalayer=1) %>%
add_overlay(place_layer)
saveRDS(scene2, "./_data/scene2.rds")
Show code
scene <- readRDS("./_data/scene.rds")
scene2 <- readRDS("./_data/scene2.rds")
scene %>% plot_3d(alt, zscale=20,
theta=30, phi=20, fov=0, zoom=0.5)
# Add polygon annotations
xyz <- sf::read_sf("~/Projects/2021-iwmi/data/mli/srtm/xyz.shp")
render_polygons(xyz, ext, data_column_top="z",
scale_data=1, color=alpha(pal["orange"], 0.8),
lit=F, light_intensity=0.01, clear_previous=T)
rglwidget()
Figure 4: Interactive 3D Scene of Selingue Dam
Threes.js
A vanilla JavaScript approach using a Three.js 3D viewer over the same area.